function Hessian = ecmnhess(Data, Covar, InvCovar, MatrixFormat)
%ECMNHESS Evaluate Hessian of negative log-likelihood function for ECMNSTD.
%	NUMPARAMS x NUMPARAMS observed Hessian matrix based on current maximum
%	likelihood parameter estimates.
%
%   Hessian = ecmnhess(Data, Covar);
%   Hessian = ecmnhess(Data, Covar, InvCovar);
%	Hessian = ecmnhess(Data, Covar, InvCovar, MatrixType);
%
% Inputs:
%   Data - NUMSAMPLES x NUMSERIES matrix of observed multivariate normal data
%		with NaNs to represent missing values.
%   Covar - NUMSAMPLES x NUMSERIES matrix of estimated covariance of Data.
%
% Optional Inputs:
%   InvCovar - Inverse of covariance matrix, i.e., inv(Covar).
%	MatrixFormat - String that identifies which parameters to be included in the
%		Hessian matrix. The default method is 'full'. The choices are:
%		'full' - (default) Compute the full Hessian matrix.
%		'meanonly' - Compute only components of the Hessian matrix associated
%			with the mean estimates.
%
% Outputs:
%   Hessian - NUMPARAMS x NUMPARAMS Hessian matrix of the observed
%		log-likelihood function based on current parameter estimates, where
%		NUMPARAMS = NUMSERIES * (NUMSERIES + 3)/2 if MatrixFormat = 'full'
%		and NUMPARAMS = NUMSERIES if MatrixFormat = 'meanonly'.
%
% WARNING: If calculating the full Hessian matrix, this routine is VERY slow.
%
% See also ECMNFISH, ECMNSTD, ECMNMLE.

%	Author(s): R.Taylor, 4-21-2005
%	Copyright 2005 The MathWorks, Inc.
%	$Revision: 1.1.6.2 $   $Date: 2005/06/17 20:23:21 $

% Step 1 - check arguments

if nargin < 2
	error('Finance:ecmnhess:MissingInputArg', ...
		'One or more of the required input arguments Data and Covar is missing.');
else
	if isempty(Data)
		error('Finance:ecmnhess:EmptyInputData', ...
			'The required input argument Data is empty.');
	end
	if isempty(Covar)
		error('Finance:ecmnhess:EmptyInputCovar', ...
			'The required input argument Covar is empty.');
	end
	
	[NumSamples, NumSeries] = size(Data);
	
	if ~all(size(Covar) == [NumSeries, NumSeries])
		error('Finance:ecmnhess:IncompatibleCovar', ...
			'The covariance matrix Covar has wrong dimensions.');
	end
end
if nargin < 3 || isempty(InvCovar) || ~all(size(InvCovar) == size(Covar))
	InvCovar = inv(Covar);
end
if nargin < 4
	MatrixFormat = 'FULL';
end

% Step 2 - initialization

MatrixFormat = upper(MatrixFormat);
if ~any(strcmp(MatrixFormat,{'MEANONLY','FULL'}))
	warning('Finance:ecmnhess:UnknownFormatString', ...
		'The MatrixFormat string is not known. Will default to FULL.');
	MatrixFormat = 'FULL';
end

if strcmp(MatrixFormat,'MEANONLY')
	NumParams = NumSeries;
else
	NumParams = NumSeries + (NumSeries * (NumSeries + 1))/2;
end

Hessian = zeros(NumParams,NumParams);

% Step 3 - main loop over data records

Map = zeros(NumSeries,1);
Count = 0;

for kk = 1:NumSamples

% Step 4 - determine and map available data in current record
    
    Map(:) = 0;
    Available = 0;
    for ii = 1:NumSeries
        if isnan(Data(kk,ii))
            Map(ii) = 0;
        else
            Map(ii) = 1;
            Available = Available + 1;
        end
    end

    if Available > 0                        % skip over empty records
        Count = Count + 1;
        
% Step 5 - construct covariance matrix subarrays

        SubCovar = Covar;
        if Available < NumSeries
            for ii = NumSeries:-1:1
                if Map(ii) == 0
                    SubCovar(:,ii) = [];
                    SubCovar(ii,:) = [];
                end
            end
            InvSubCovar = inv(SubCovar);
        else
            InvSubCovar = InvCovar;
        end
        
% Step 6 - do partials wrt Mean for current data record

        ii = 0;
        for i = 1:NumSeries
            if Map(i) > 0
                ii = ii + 1;
            end
            
            jj = 0;
            for j = 1:i;
                if Map(j) > 0
                    jj = jj + 1;
                end
                
                if (Map(i) > 0) && (Map(j) > 0)
                    Hessian(i,j) = Hessian(i,j) + InvSubCovar(ii,jj);
                    Hessian(j,i) = Hessian(i,j);
                end
            end
        end

% Step 7 - do partials wrt Covar for current data record

		if strcmp(MatrixFormat,'FULL')

			GradC1 = zeros(Available,Available);
			GradC2 = zeros(Available,Available);

			p1 = 0;
			i = NumSeries;
			for i1 = 1:NumSeries
				if Map(i1) > 0
					p1 = p1 + 1;
				end

				q1 = 0;
				for j1 = 1:i1
					i = i + 1;

					if Map(j1) > 0
						q1 = q1 + 1;
					end

					if (Map(i1) > 0) && (Map(j1) > 0)
						GradC1(p1,q1) = 1.0;
						GradC1(q1,p1) = 1.0;
					end

					p2 = 0;
					j = NumSeries;
					for i2 = 1:NumSeries
						if Map(i2) > 0
							p2 = p2 + 1;
						end

						q2 = 0;
						for j2 = 1:i2
							j = j + 1;

							if Map(j2) > 0
								q2 = q2 + 1;
							end

							% dC/dtheta(i) = dC/dC(i1,j1)
							% dC/dtheta(j) = dC/dC(i2,j2)

							if (j <= i) && (Map(i1) > 0) && (Map(j1) > 0) && (Map(i2) > 0) && (Map(j2) > 0)
								GradC2(p2,q2) = 1.0;
								GradC2(q2,p2) = 1.0;

								Temp1 = InvSubCovar*GradC1;
								Temp2 = InvSubCovar*GradC2;

								Hessian(i,j) = Hessian(i,j) + 0.5*trace(Temp1*Temp2);
								Hessian(j,i) = Hessian(i,j);

								GradC2(p2,q2) = 0.0;
								GradC2(q2,p2) = 0.0;
							end
						end
					end

					if (Map(i1) > 0) && (Map(j1) > 0)   % undo dC/dtheta(i)
						GradC1(p1,q1) = 0.0;
						GradC1(q1,p1) = 0.0;
					end
				end
            end
        end
    end
end

% Step 8 - normalize hessian

Hessian = Hessian ./ Count;
